This notebook is a simple notebook to show the loading, plotting on a map and gridding Darwin CPOL data to a Cartesian grid at two different resolutions
In [27]:
#load relevant modules, set plotting to inline
import pyart
from matplotlib import pyplot as plt
import numpy as np
import wradlib
from pyart.core.echo_classes import *
from matplotlib import rc
from netCDF4 import num2date
print pyart.__version__
%matplotlib inline
In [48]:
#Data is at: https://engineering.arm.gov/~collis/cpol/
In [3]:
filename = '/data/cpol/Gunn_Pt_rapic_20120315_102002_PPI.UF'
radar = pyart.io.read(filename)
In [5]:
print radar.fields.keys()
In [20]:
myd = pyart.graph.RadarMapDisplay(radar)
fig = plt.figure(figsize = [18,10])
myd.plot_ppi_map('corrected_reflectivity', vmin = -16, vmax = 64, resolution = 'h')
m= myd.basemap
runway1x,runway1y=m(np.array([130.865321, 130.89427]),np.array([-12.40945, -12.419184]))
runway2x,runway2y=m(np.array([130.870643, 130.870696]), np.array([-12.409041, -12.422505]))
apannotx, apannoty=m(130.75, -12.37)
armx, army=m(130.891848,-12.424575)
cpolx, cpoly=m(131.044402,-12.249027)
bayx, bayy=m(130.9,-12.3)
m.drawcoastlines()
m.drawmapboundary(fill_color='blue')
m.drawparallels(np.linspace(-11, -13, 11),labels=[1,0,0,0])
m.drawmeridians(np.linspace(130,132, 11),labels=[0,0,0,1])
m.drawrivers()
m.plot(runway1x, runway1y, 'k-')
m.plot(runway2x, runway2y, 'k-')
# P.arrow( x, y, dx, dy, **kwargs )
#P.arrow( 0.5, 0.8, 0.0, -0.2, fc="k", ec="k",
plt.arrow(apannotx, apannoty, runway1x[0]-apannotx-4000,
runway1y[0]-apannoty+1000, head_width=3000,
head_length=3000, fc='w', ec='w' )
plt.text(apannotx-5000,apannoty,'Airport', fontdict={'color':'white', 'weight':'bold'})
m.plot(armx,army,'bo')
plt.text(armx+1500, army-3000, 'ARM Darwin\n site', fontdict={'color':'black', 'weight':'bold'})
m.plot(cpolx,cpoly,'bo')
plt.text(cpolx+1500, cpoly+1000, 'CPOL Radar', fontdict={'color':'black', 'weight':'bold'})
plt.text(bayx,bayy,'Shoal \n Bay', fontdict={'color':'white', 'weight':'bold'})
m.drawmapscale(130.8, -13.2, 131, -13.2, 100, barstyle='fancy', fontcolor='w', fillcolor1='w', fillcolor2='k')
Out[20]:
In [21]:
cpol_map = pyart.map.grid_from_radars((radar,),
grid_shape=(35, 401, 401),
grid_limits=((0, 17000), (-150000, 150000), (-150000, 150000)),
grid_origin = (-12.249027,131.044402),
fields=radar.fields.keys(),
refl_field='corrected_reflectivity',
max_refl=100.)
In [22]:
pyart.io.write_grid('/data/cpol/cpol_20120315_102002.nc', cpol_map)
In [47]:
display = pyart.graph.GridMapDisplay(cpol_map)
# create the figure
font = {'size': 14}
rc('font', **font)
fig = plt.figure(figsize=[15, 8])
# panel sizes
map_panel_axes = [0.05, 0.05, .4, .80]
x_cut_panel_axes = [0.55, 0.10, .4, .30]
y_cut_panel_axes = [0.55, 0.50, .4, .30]
colorbar_panel_axes = [0.05, 0.90, .4, .03]
# parameters
level = 5
vmin = -16
vmax = 64
lat = -11.89
lon = 130.7
# panel 1, basemap, radar reflectivity a
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution = 'h')
display.plot_grid('corrected_reflectivity', level=level, vmin=vmin, vmax=vmax)
display.plot_crosshairs(lon=lon, lat=lat)
m= display.basemap
runway1x,runway1y=m(np.array([130.865321, 130.89427]),np.array([-12.40945, -12.419184]))
runway2x,runway2y=m(np.array([130.870643, 130.870696]), np.array([-12.409041, -12.422505]))
apannotx, apannoty=m(130.75, -12.37)
armx, army=m(130.891848,-12.424575)
cpolx, cpoly=m(131.044402,-12.249027)
bayx, bayy=m(130.9,-12.3)
m.drawcoastlines()
m.drawmapboundary()
m.drawparallels(np.linspace(-11, -13, 11),labels=[1,0,0,0])
m.drawmeridians(np.linspace(130,132, 6),labels=[0,0,0,1])
m.drawrivers()
m.plot(runway1x, runway1y, 'k-')
m.plot(runway2x, runway2y, 'k-')
# P.arrow( x, y, dx, dy, **kwargs )
#P.arrow( 0.5, 0.8, 0.0, -0.2, fc="k", ec="k",
plt.arrow(apannotx, apannoty, runway1x[0]-apannotx-4000,
runway1y[0]-apannoty+1000, head_width=3000,
head_length=3000, fc='w', ec='w' )
plt.text(apannotx-5000,apannoty,'Airport', fontdict={'color':'white', 'weight':'bold','size': 12})
m.plot(armx,army,'bo')
plt.text(armx+1500, army-3000, 'ARM Darwin\n site', fontdict={'color':'black', 'weight':'bold','size': 12})
m.plot(cpolx,cpoly,'bo')
plt.text(cpolx+1500, cpoly+1000, 'CPOL Radar', fontdict={'color':'black', 'weight':'bold', 'size': 12})
plt.text(bayx,bayy,'Shoal \n Bay', fontdict={'color':'white', 'weight':'bold', 'size': 12})
m.drawmapscale(130.8, -13.0, 131, -13.0, 100, barstyle='fancy', fontcolor='b', fillcolor1='b', fillcolor2='k')
# colorbar
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(cax=cbax)
# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from SGP CF (km)')
# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)
# add a title
slc_height = cpol_map.axes['z_disp']['data'][level]
dts = num2date(cpol_map.axes['time']['data'], cpol_map.axes['time']['units'])
datestr = dts[0].strftime('%H:%M Z on %Y-%m-%d')
title = 'Sliced at ' + str(slc_height) + ' meters at ' + datestr
fig.text(0.5, 0.9, title)
Out[47]:
In [41]:
cpol_map_hr = pyart.map.grid_from_radars((radar,),
grid_shape=(35, 601, 601),
grid_limits=((0, 17000), (-150000, 150000), (-150000, 150000)),
grid_origin = (-12.249027,131.044402),
fields=radar.fields.keys(),
refl_field='corrected_reflectivity',
max_refl=100.)
In [44]:
pyart.io.write_grid('/data/cpol/cpol_hr_20120315_102002.nc', cpol_map_hr)
In [46]:
display = pyart.graph.GridMapDisplay(cpol_map_hr)
# create the figure
font = {'size': 14}
rc('font', **font)
fig = plt.figure(figsize=[15, 8])
# panel sizes
map_panel_axes = [0.05, 0.05, .4, .80]
x_cut_panel_axes = [0.55, 0.10, .4, .30]
y_cut_panel_axes = [0.55, 0.50, .4, .30]
colorbar_panel_axes = [0.05, 0.90, .4, .03]
# parameters
level = 5
vmin = -16
vmax = 64
lat = -11.89
lon = 130.7
# panel 1, basemap, radar reflectivity a
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution = 'h')
display.plot_grid('corrected_reflectivity', level=level, vmin=vmin, vmax=vmax)
display.plot_crosshairs(lon=lon, lat=lat)
m= display.basemap
runway1x,runway1y=m(np.array([130.865321, 130.89427]),np.array([-12.40945, -12.419184]))
runway2x,runway2y=m(np.array([130.870643, 130.870696]), np.array([-12.409041, -12.422505]))
apannotx, apannoty=m(130.75, -12.37)
armx, army=m(130.891848,-12.424575)
cpolx, cpoly=m(131.044402,-12.249027)
bayx, bayy=m(130.9,-12.3)
m.drawcoastlines()
m.drawmapboundary()
m.drawparallels(np.linspace(-11, -13, 11),labels=[1,0,0,0])
m.drawmeridians(np.linspace(130,132, 6),labels=[0,0,0,1])
m.drawrivers()
m.plot(runway1x, runway1y, 'k-')
m.plot(runway2x, runway2y, 'k-')
# P.arrow( x, y, dx, dy, **kwargs )
#P.arrow( 0.5, 0.8, 0.0, -0.2, fc="k", ec="k",
plt.arrow(apannotx, apannoty, runway1x[0]-apannotx-4000,
runway1y[0]-apannoty+1000, head_width=3000,
head_length=3000, fc='w', ec='w' )
plt.text(apannotx-5000,apannoty,'Airport', fontdict={'color':'white', 'weight':'bold','size': 12})
m.plot(armx,army,'bo')
plt.text(armx+1500, army-3000, 'ARM Darwin\n site', fontdict={'color':'black', 'weight':'bold','size': 12})
m.plot(cpolx,cpoly,'bo')
plt.text(cpolx+1500, cpoly+1000, 'CPOL Radar', fontdict={'color':'black', 'weight':'bold', 'size': 12})
plt.text(bayx,bayy,'Shoal \n Bay', fontdict={'color':'white', 'weight':'bold', 'size': 12})
m.drawmapscale(130.8, -13.0, 131, -13.0, 100, barstyle='fancy', fontcolor='b', fillcolor1='b', fillcolor2='k')
# colorbar
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(cax=cbax)
# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from SGP CF (km)')
# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)
# add a title
slc_height = cpol_map_hr.axes['z_disp']['data'][level]
dts = num2date(cpol_map_hr.axes['time']['data'], cpol_map_hr.axes['time']['units'])
datestr = dts[0].strftime('%H:%M Z on %Y-%m-%d')
title = 'Sliced at ' + str(slc_height) + ' meters at ' + datestr
fig.text(0.5, 0.9, title)
Out[46]:
In [ ]: